Experimental Quantum Simulation of Entanglement in Many-body Systems 
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We employ a nuclear magnetic resonance (NMR) quantum information processor to simulate 
the ground state of an XXZ spin chain and measure its NMR analog of entanglement, or pseudo- 
entanglement. The observed pseudo-entanglement for a small-size system already displays singular- 
ity, a signature which is qualitatively similar to that in the thermodynamical limit across quantum 
phase transitions, including an infinite-order critical point. The experimental results illustrate a 
successful approach to investigate quantum correlations in many-body systems using quantum sim- 
ulators. 
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Entanglement, delineated as the non-local correlation, 
is one "spooky" characteristic trait of quantum mechan- 
ics The famous dispute between Bohr and Ein- 
stein on the fundamental question of quantum mechanics, 
Schrodingcr-cat paradox, and the transitions from quan- 
tum to classical worlds essentially involve entanglement. 
Recent development of quantum information rekindles 
the interest in entanglement, more as a possible resource 
for information processing j2|. Various methods have 
been proposed to characterize entanglement qualitatively 
and quantitatively Q . One immediate application of the 
entanglement is the investigation of quantum phase tran- 
sitions (QPTs) li-Q in many-body systems, which oc- 
cur at zero temperature (T=0 K), where the transitions 
are driven by quantum fluctuations and the ground-state 
wave function is expected to develop drastic change. The 
entanglement properties extend and complement the tra- 
ditional statistical-physical methods for QPTs, such as 
the correlation functions and low-lying excitation spec- 
tra. However, how to describe and measure entanglement 
in many-body systems is still a challenging task in both 
theoretical and experimental aspects Most schemes 
for directly measuring entanglement focus on the entan- 
glement between two qubits [3, • Although the degree 
of entanglement for a medium-size or lager system can in 
principle be probed it has not been experimentally 
measured directly. 

In contrast to classical approaches, quantum simula- 
tors [l^ provide a promising approach for investigating 
many-body systems and enable one to efficiently simu- 
late other quantum systems by actively controlling and 
manipulating a certain quantum system, and to test, 
probe and unveil new physical phenomena. One inter- 
esting aspect is to simulate the ground states of many- 
body systems, where usually rich phases can exist, such 
as ferromagnetism, superfluidity, and quantum Hall ef- 
fect, just to name a few. In this Letter we experimentally 
simulate the ground state of an XXZ spin chain [l3| in 
a liquid-state nuclear magnetic resonance (NMR) quan- 
tum information processor fl^ and directly measure a 



global multipartite entanglement — the geometric entan- 
glement (GE) [l^, [l^ in a version of NMR analog, or 
pseudo-entanglement. Exploiting the probed behavior 
of GE, we identify two QPTs, which in the thermody- 
namic limit correspond to the first- and oo-ordcrs, respec- 
tively. In the (X)-order QPT, also known as Kosterlitz- 
Thouless (KT) transition HI] , the ground-state en- 
ergy is not singular. Consequently the detection of the 
critical point in the KT transition niaypose a challenge 
for the correlation-based approaches Ia,[i9| , which rely on 
the singularity of ground-state energy. Surprisingly, the 
GE turns out to be non-analytical but of different types 
of singularity at the first- and 00-order transitions [20|. 
Remarkably, the qualitative features of the ground-state 
entanglement in the thermodynamic limit displayed near 
both transitions persist even for a small-size system, on 
which our experiment is performed. 

The GE of a pure many-spin quantum state 1$) 
is captured by the maximal overlap [l^, [l^ A^ax = 
max* K^fl^)!, and is defined as i^iog^ = — logjA^^^, 

where l^*) = ^f^i IV"^*^) denotes all product (i.e. unen- 
tangled) states of the iV-spin system. From the point of 
view of local measurements, the GE is essentially (mod- 
ulo a logarithmic function) the maximal probability that 
can be achieved by a local projective measurement on 
every site, and the closest product state signifies the op- 
timal measurement setting. The GE has been employed 
to study QPTs [l^ [2l|, local state discrimination |22j 
and entanglement as computational resources [23| . 
The XXZ spin chain is described by the Hamiltonian 
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Hxxz = Y.{X,X^+l + Y,Y,+, + jZ,Z,+,) (1) 
1=1 

where A,;, Yi, Zi denote the Pauli matrices with i indi- 
cating the spin location, and 7 is the control parameter 
for QPTs. We use the periodic boundary condition with 
A^ -I- 1 = 1. The XXZ chain can be exactly solved by the 
so-called Bethe Ansatz and exhibits rich phase diagrams 
in the ground state [l3| . In the thermodynamic limit. 
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FIG. 1: (Color online). Theoretical results in the 4-spin XXZ 
chain, (a) Energy level of the ground state, (b) The overlap 
square A^(7) and entanglement £"105,(7). The solid, dashed 
and dash-dotted curves show Af (7) for 1*1) = |1111), 1*2) = 

I H \ — ) and {'^3) = |0101), respectively. £1052(7) is shown 

as the dotted curve. The jump at 7 = — 1 and the cusp at 
7 = 1 in -Biog2 (7) indicate the transition points for the QPTs, 
with the first- and 00- orders, respectively. 



glcment features pertinent to the QPTs in the thermo- 
dynamic limit will survive. The ground-state energy and 
wave function of the 4-spin chain are represented as (see 
Supplemental Material) 
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where a € (— 7r/2, 0) is given via tan(2a) = —2\/2/j, and 
= (10101) + |1010))/V2, \ct>2) = (IllOO) + 10011) + 
1 1001) + |0110))/2. Fig. [U (a) shows Eg as a function 
of 7. One should note that Eg is continuous at 7 = — 1 
while the ground state \g) is discontinuous. 

In order to obtain the ground-state GE, one needs to 
search for the closest product state j^**) In fact, 

we can choose the product states 
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to obtain the corresponding entanglement in the respec- 
tive range of 7. To anticipate the experimental proce- 
dure, we shall measure the ground-state overlap listed as 



for 7 < —1, the system has the ground state with the 
ferromagnetic (FM) Ising phase. At 7 = —1 a first-order 
QPT occurs. For — 1 < 7 < 1, the system is in a gapless 
phase or XY-like phase. At 7 = 1 there is an 00-order 
or a KT transition (T^j, where the ground-state energy, 
however, is analytic across the transition, and so is any 
correlation function. For 7 > 1, the system is in the Neel- 
like antiferromagnetic (AFM) phase. The ground state 
is asymptotically doubly degenerate. However, the exci- 
tations above the ground space have a gap. For 7^1, 
the ground state takes a Nccl or Ising AFM state i.e., 
|... 10101010...), where |0) = | t) and |1) = | ;). 

In the XXZ chain ([T]), the GE displays a jump across 
7 = — 1 but a cusp (i.e., the derivative is discontinuous) 
across 7 = 1 [2^. Both features are present for small- 
size systems, as well as in the thermodynamic limit. For 
7 < — 1 the GE is essentially zero. Regardless of the sys- 
tem size (as long as it is even), for —1 < 7 < 1, the closest 

product state is found to be | H \ — ...), whereas for 

7 > 1, the closest product state is found to be |0101...), 
where |±) = (|0) ± \l))/V2 Right at the KT point 
7 = 1, because of rotational symmetry, the closest prod- 
uct states are |0 (/)"*" ^(/)^...), where \(j>) and are any 
arbitrary orthonormal qubit states. The singular behav- 
ior of ground-state GE can be used to probe the KT 
transition, and there is no need to know the low-lying 
spectrum. 

In implementation we use a 4-spin chain. The entan- 



Ai(7) = (*i|,9)= ' , (5) 



A-2h) = {^-M = 1 . JJ (6) 

I cosa — ^ sma (7 > — 1) 

fo (7<-l) 
^^(^)-<*^l^>-|^cosa (^>-l)- 

From Eqs. ^ and ([7]), one finds that A2(7) and A3 (7) 
cross at 7 = 1. Fig. [T](b) shows the theoretical prediction 
for A^(7) {i = 1, 2, 3) and the entanglement Ei^g^. The 
jump in the entanglement at 7 = — 1 and the cusp at 
7 = 1 signify the two QPT points. 

In exp eriment, we choose the four carbons in crotonic 
acid [23] dissolved in d6-acetone as the four qubits. We 
generate the ground states using quantum networks and 
implement the quantum gates by GRAPE pulses [26| . 
Various ground-states can be created by varying the sin- 
gle spin rotations that can be easily implemented (see 
Supplemental Material). In principle, one can employ 
an iterative method to experimentally measure the GE 
of the ground state (see Supplemental Material). To 
demonstrate the proof-of-principle simulation of quan- 
tum entanglement, instead, we first measure the overlap 
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FIG. 2: (Color online). Experimentally measured Af(7) (a) 
and -Eiog2 (b) for various 7. In figure (a), the experimen- 
tal data are shown as "x" and for A- = K*^|5)|^ 
corresponding to I*!*!), |*I'2) and l^a). The measured A^iij) 
for 7 < — 1 are indicated as the dotted, dashed and dash- 
dotted lines. In the region 7 > —1, Af (7) can be fitted as the 
dotted line. Through fitting the points for A2(7) and A|(7) 
using polynomial functions, we obtain the two solid curves 
that cross at point 7 — 0.92, close to the theoretical point at 
7=1. The thick dashed and dash-dotted curves show the 
fitting results using the theoretical A2(7) and A|(7) by in- 
troducing decay factors 0.69 and 0.71, respectively. In figure 
(b), in range 7 < — 1, Eio^^ is shown as "-|-". For 7 > — 1, 
we rescale the measured A2(7) and A3 (7) as A2(7)/0.69 and 
A3 (7)70.71, respectively. The dotted and dashed curves that 
cross at 7 = 1.02 show the fitting results of the rescaled data 
using polynomial functions. The expected -Eiogj after rescal- 
ing is indicated by "o". 



of the ground state with several product states ([4]), vifhich 
contain the closest product states. From the measure- 
ment with the already known closest product states, we 
can obtain the ground-state GE. Next, to show that the 
obtained results are the optimum, we vary the product 
states to test the optimality. 

The experimentally measured Af (7) for various 7 are 
shown in Fig. [2] (a). The measured A^(7) for 7 < — 1 are 
0.92, 0.048 and 0.019, indicated as the dotted, dashed and 
dash-dotted lines. The corresponding theoretical values 
are 1, 1/16 and 0, respectively. In the region 7 > — 1, 
Af (7) can be fitted as A^(7) = 0.014, shown as the dotted 
line, corresponding to in theory. 

We perform polynomial fits to the measured A^ (7) and 
A|(7), and obtain the two solid curves that cross at the 
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FIG. 3: (Color online). Theoretical (a) and experimen- 
tally measured (b) A^ for various product states created as 
C/p(^)|0101) [see Eq. (gj] . Three ground states for 7 = -0.9, 
1 and 3 are chosen and the corresponding A'^ are shown as 
the the solid, dash-dotted and dashed curves in figure (a), 
respectively. The experimental data are shown as "■" 
and "x" for 7 — —0.9, 1 and 3 in figure (b), respectively. In 
comparison with the theoretical values, they can be fitted as 
0.66A^ 0.68A^ and 0.71A^ shown as the solid, dash-dotted 
and dashed curves. 



point 7 = 0.92, which is very close to the theoretically 
predicted transition point at 7 = 1. The discrepancy be- 
tween experiment and theory mainly comes from the dif- 
ferent experimental errors in measuring A2(7) and A§(7). 
The jump at 7 = — 1 and the cusp at 7 = 0.92 reflect the 
different types of QPT points. 

In order to faithfully estimate the performance of 
the experiment in measuring A|(7) and A|(7) in the 
range 7 > —1, we introduce two decay factors 02 
and aa to fit the experimental data as [A2 3(7)]oxp = 
Q!2,3[A| 3(7)]thoory, shown as the thick dashed and dash- 
dotted curves in Fig. [2] (a) with the best scale-factors 
as a2 = 0.69 and as = 0.71, respectively. The differ- 
ence between the decay factors comes from the different 
operations in measuring A2(7) and A3 (7). In Fig. [2] 
(b), we exploit the decay factors to rescale experimen- 
tal values of [A2 3]oxp/q;2,3, from which we obtain the ex- 
pected values of pseudo-entanglement shown as " o" . The 
rescaled — log2([A2 3]gxp/Q!2.3) can be fitted as the dotted 
and dashed curves that cross at 7 = 1.02. 

In principle, we do not need to know the closest prod- 
uct states in order to measure the entanglement. In the 
Supplemental Material, we describe an iterative proce- 
dure to search for them and this procedure can be im- 
plemented in experiment. For proof-of-principlc demon- 
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stration of the optimality experimentally, we simplify the 
procedure and vary the product states |^'(/3)) by 

=C/p(/?)|0101), (8) 

where Up{(3) = (S'j=i e~*^^^/^, and experimentally mea- 
sure = |(^'(/3)|(7)p for various /3 at three different lo- 
cations of the phase diagram, corresponding to 7 = —0.9, 
1. and 3, respectively. The theoretical and experimental 
results arc shown in Figs.[3lja) and (b), respectively. The 
experimental data are compared to the theoretical values 
of 0.66A2, 0.68A2 and O.TlA^, shown in Fig.[3](b). One 
finds that the maximum of A^ occurs at /3 = 7r/2 and 
for 7 = —0.9 and 3, respectively. These correspond to 

the respective closest product states, {'^2) = I H i — ) 

and I'I's) = |0101), predicted theoretically. Remarkably, 
for 7 = 1, where the cxj-order QPT occurs, A^ is a con- 
stant independent of /3, as we have expected and noted 
earlier. This also means that arbitrary states prepared 
by Eq. ^ can be chosen to measure the entanglement 
at 7 = 1, and this gives additional confirmation that the 
created state at the KT point is rotationally invariant. 

The experiment duration of the preparation of the 
ground states for 7 > —1 is about 160 ms, which is 
non-negligible (about 17%) compared to the coherence 
time T2. Consequently the decay of the signals due to 



the limitation of coherence time is one of main sources of 
errors. Additionally the imperfection of pulses and inho- 
mogeneities of magnetic fields also contribute to errors. 
The deviations of the experimental data from the the- 
oretical fitting in Fig. [3] (b) represent the effects of the 
errors that depend on the rotation angles, or the product 
states. In particular the fluctuation of the data for 7 = 1 
in Fig. [3] (b) confirms the explanation for the shift of the 
measured cusp in Fig. [2] (a) . 

In conclusion we demonstrate the non-analytic proper- 
ties of many-body systems in a quantum simulator using 
NMR. The QPTs with first- and 00-orders in the XXZ 
spin chain are detected by directly measuring the pseudo- 
entanglement of the ground states created by quantum 
gates. An alternative approach for creating ground states 
would be via adiabatic evolution [27j . Our prelimi- 
nary numerical analysis indicates that ground states for 
7 > — 1 can be approximately generated with high fi- 
delity (e.g. > 0.998) by the adiabatic evolution from the 
ground state at a large 7. The experimental implemen- 
tation is a possible future direction. 
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Supplemental IMaterial 

Appendix A: Method for computing and measuring 
the maximal overlap by iteration 

Here we describe an iterative method to compute 
the maximal overlap, of which the motivation comes 
from the density-matrix-renormalization-group (DMRG) 
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or matrix-product-state (MPS) variational method 
This method can not only be implemented numerically, 
but can also be carried out experimentally. To compute 
the maximal overlap for the state with respect to 
product states l^*) = 0^=1 we use the Lagrange 

multiplier A to enforce the constraint (^'1^') = 1, 

/(*)^(*|g)(5|*)-A(*|*). (Al) 

Maximizing / with respect to the local product state 
l^/!*^*)), we obtain the extremal condition 

H^^^IV^'')) =A7VW|^W), (A2) 

where n^!^ = (0f^,(V^(^')|)|.9>(5l((8)f^. 1^^^')) is propor- 
tional to a local projector (labeled by 1 1/)^*^ | at i- 

sitc), and the normalization A^''^ = (^^^^(V'^"''' IV''"'"'): is 
unity if all the local states are properly normalized. From 
the viewpoint of the variational MPS, one fixes all local 
states IV'^-'') but IV'^''') and solves for the corresponding 
optimal IV'*'*'') and repeats the same procedure for i + 1, 
1 + 2, etc. until the A^-th site and sweeps the procedure 
back and forth until the eigenvalue A converges. The 
converged value |Ap is the square of the maximal over- 
lap A^^,,. 

Experimentally, this procedure means that one chooses 
randomly the local measurement basis and picks arbitrar- 
ily one direction (i.e., rank-one projector) for each site, 
say, IV'*'"''') for the j-th site and only varies the basis for 
one site, say, i-th at a time with all others fixed until one 
reaches a basis where the measurement outcome along 
one direction occurs with the most probability. Then one 
moves to the next site, say, (i + l)-th site, and finds the 
optimal direction and repeats this procedure by sweeping 
back and forth until the probability for the most likely 
outcome converges. 

Numerically, it appears that one has to solve the above 
eigenvalue problem. But as the effective local Hamilto- 
nian is a projector, one immediately sees that the optimal 
rank-one projector is exactly the one ( | |) given 
by H^J^. One thus replaces by and repeats the 
procedure at other sites until the overlap converges. Ex- 
perimentally, the apparatus setting in line with the pro- 
jector gives the local maximum output probability. The 
search for the optimal direction at the «-th site need not 
be a blind search, as a tomography (conditioned on all 
other sites being measured in their respective |'0^'''')) will 
enable the determination of the optimal local direction 
10^')). The whole procedure, either numerically or ex- 
perimentally, thus becomes an iterative procedure, given 
an initial choice of {jV'^-'-')}. We have performed such a 
numerical procedure and have found that this procedure 
converges efficiently to the maximal overlap. The conver- 
gence for the ground state of four-spin chain is very fast 
and the result is accurate; see Fig. |31 In our experiment, 
we implement a simplified version to directly measure the 
ground-state entanglement. 
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10 20 30 40 50 
Iteration step 

FIG. 4: (Color online). Numerical implementation of the 
iteration procedure for the maximal overlap in the 4-spin XXZ 
chain. The overlap square is shown as marked as "•" and 
"x" for 7 = —0.9, 1 and 3, respectively. An initial 4-spin 
product state is randomly chosen. One round of sweep starts 
from the first spin, reaches the forth spin, and sweeps back 
to the second spin. One round thus contains six steps. Eight 
iteration rounds are shown. The maximal overlap squares 
converge rapidly to their exact values, which are indicated by 
the thin lines, after just a few rounds. We have checked that 

the converged product states are | H 1 — ) (up to a relative 

phase), (with {(f>'''\(j)) = 0), and |0101), respectively. 

Appendix B: Method for solving the 4-spin XXZ 
chain 

For 7 < — 1 the ground state is doubly degenerate 
with ground states being |0000) and |1111). To avoid 
the complication due to the degeneracy, we can intro- 
duce an additional small Zeeman term = 
with < -Bz ^ 1 into the Hamiltonian to lift the energy 
of |0000) so that the ground state \g) becomes |1111). 
The small Zeeman term does not affect the universality 
class of phase transition of the ground state, because it 
commutes with Hxxz- 

For 7 > — 1, the ground state is not a simple product 
state. Due to the absence of an external field, the conser- 
vation of z-component total angular momentum and the 
periodic translation invariance give rise to only two rele- 
vant "Bcthe-ansatz" basis states for the ground states: 

|0i) ^ -L(|oiOI) + 11010)) (Bl) 

102) = i(|1100) + 10011) + 11001) + 10110)). (B2) 

By solving the effective Hamiltonian in the subspace 
spanned by {|(/)i), |02)} 
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we obtain that the ground state is 

\g) ~ \(f>i) cosa + 102) sinck (B4) 
where a € (— 7r/2, 0) is given via 

tan(2a) -2\/2/7, (B5) 
and that the ground energy is 

(7) = -27-2^7^ + 8. (B6) 

Appendix C: Method for experimental 
implementation 

The experiment is performed in a Bruker DRX 700 
MHz spectrometer. The structure of the molecule of cro- 
tonic acid and the parameters of the four spin qubits are 
shown in Fig. [5] (a). The protons are decoupled in the 
whole experiment. The initial pseudo-pure state |0000) 
is prepared by spatial averaging 0, Q , and chosen as the 
reference state for normalizing the signals in the following 
partial state tomography. 

The quantum circuit shown as Fig. [5] (b) illustrates 
the experiment protocol for 7 > —1. The ground state 
is created by U3U2{a)Ui, indicated by the three dashed 
blocks, respectively. We optimize Ui and C/3, which are 
independent of a, as two long (40 ms duration) GRAPE 
pulses, respectively, where GRAPE stands for gradient 
ascent pulse engineering. The theoretical fidelity for Ui 
and C/3 is larger than 99%. To save time in searching 
GRAPE pulses, we further decompose 1/2(0) into simple 
gates shown as the sequence in Fig. [5] (c), where each 
gate is implemented by one GRAPE pulse. The duration 
of the pulse for the spin coupling evolution is 20 ms, and 
the duration of other pulses is 0.5 ms. The theoretical 
fidelity for each pulse in Fig. [5] (c) is larger than 99.5%. 



An arbitrary ground state for 7 > — 1 can be generated 
through varying a in the single spin operation, which is 
much easier to find than U2 (a) in the GRAPE algorithm. 
For the case of 7 < —1, we replace C/3C/2(a)C/i by four 
NOT gates implemented by four tt pulses applied to the 
four qubits respectively to create the ground state |1111) 
from 10000). 

To measure the overlap between \g) and an arbitrary 
product state |^), we re- write the overlap A = (^'l^) in 
form of 

A = {b\Ul\g) (CI) 

where \b) denotes a computational basis and Up\b) = l^*) 
3 . Here we choose Up as 

4 

f/p(/?) =(g)e-''3^^/2. (C2) 

Since \'i>i) and [^'3) are already the computational basis, 
we can simply choose \b) as |^'i) and |5'3), respectively, 
and take Up as the identity operation by setting /3 = 0, 
for obtaining Ai and A3 from Eq. (|Cip . |'I'2) is not a 
computational basis. We can choose \b) = |0101) and /3 = 
tt/2 for obtaining A2, noting that 1^-2) = [/p(7r/2)|0101). 
In the density-matrix form, Eq. (|C1[) is represented as 

A'=Tri\b){b\p) (C3) 

where p = U^{\g){g\)Up. From Eq. (|C3|) . one finds that 
A^ is encoded as the diagonal element \b){b\ of p. We 
exploit phase cycling to remove all the non-diagonal ele- 
ments, and then reconstruct all the diagonal terms of p 
using partial state tomographgy through four tt/2 read- 
out pulses selective for C1-C4, respectively. A^ is there- 
fore extracted from the diagonal terms. 
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FIG. 5: (Color online). Experimental protocol, (a) The 
molecule structure (inset) and the parameters of the four 
carbon-13 spins. The diagonal terms show the chemical shifts 
and the non-diagonal terms show the strength of J- couplings 
in Hz. The transversal relaxation times T2 measured by a 
Hahn echo are listed in the rightmost column, (b) Quan- 
tum circuit for creating the ground state and measuring its 
overlap A; (7) for 7 > —1. Here X denotes a NOT gate. 
R±a = exp(±iQy), and H ( H*) denotes a gate transform- 
ing state |0) (|1)) to (|0) + \l))/V2, and |1) (|0)) to (|0) - 
\1))/V2, respectively. The other gates are illustrated in the 
legend box. U2 is decomposed as the gate sequence shown in 
(c), for creating an arbitrary ground state for 7 > — 1 through 
varying a [see Eq. ()B5|) ]. Up denotes a rotation for obtaining 
the overlap of the ground-state with Up\b), where \b) denotes 
a computational basis. In dephasing operation, j, k, m, n 
are chosen as and 1, respectively, to average out all non- 
diagonal terms in the density matrix to zero. Tomograhpy of 
the diagonal terms requires four 7r/2 pulses selective for CI 
to C4, respectively. 



